Elizabethkingia (oligotyping)

Load packages, paths, functions

# Load main packages, paths and custom functions
source("../../../source/main_packages.R")
source("../../../source/functions.R")

# Load supplementary packages
packages <- c("RColorBrewer", "ggpubr", "cowplot", "Biostrings", "openxlsx", "kableExtra")
invisible(lapply(packages, require, character.only = TRUE))

Preparation

Tables preparation

Seqtab

# move to oligotyping directory
path_elizabeth <- "../../../../output/2_Oligotyping/2A/Elizabethkingia/2A_oligotyping_Elizabethkingia_sequences-c1-s1-a0.0-A0-M10"

# load the matrix count table
matrix_count <- read.table(paste0(path_elizabeth, "/MATRIX-COUNT.txt"), header = TRUE) %>% t()

# arrange it
colnames(matrix_count) <- matrix_count[1,]
matrix_count <- matrix_count[-1,]
matrix_count <- matrix_count %>% as.data.frame()

# print it
matrix_count %>%
  kbl() %>%
  kable_paper("hover", full_width = F)

	CTC1	CTC10	CTC11	CTC12	CTC13	CTC14	CTC15	CTC2	CTC3	CTC4	CTC5	CTC6	CTC7	CTC9	NP27	NP30	NP36	S126	S146	S147	S154	S160	S163	S164	S165	S166	S18	S19	S20	S21	S23	S24	S25	S26	S27	S28	S29	S30	S31	S32	S33	S34	S35	S36	S37	S38	S39	S40	S41	S44	S45	S46	S47	S48	S49	S50	S51	S52
C	19626	592	5363	287	2687	1000	14770	1700	2520	1075	326	10704	15	2688	11	1	6	1	4	1	1	1	1	6	2	2	1	12104	12	6	33	53	3	8	1755	5	4	1933	4229	2020	2060	225	392	17944	3328	4744	4931	4388	1	36	32	9	6	221	43	71	61	47
G	3017	80	963	43	437	150	2088	254	378	165	60	1649	5	399	1	0	1	0	0	0	0	0	0	1	0	4	0	277	0	0	5	1	0	1	27	0	1	55	121	36	66	5	10	390	71	73	110	130	0	1	1	0	1	6	2	3	1	3

Taxonomy

# load the fasta table
fasta <- readDNAStringSet(paste0(path_elizabeth, "/OLIGO-REPRESENTATIVES.fasta"))

# arrange it
fasta <- fasta %>% as.data.frame()
colnames(fasta) <- "seq"
fasta$oligotype <- rownames(fasta)
fasta <- fasta %>% dplyr::select(-c(seq))

# print it
fasta %>%
  kbl() %>%
  kable_paper("hover", full_width = F)

	oligotype
C	C
G	G

Change oligotype name by oligotype / MED nodes in the matrix count

# Reference file 

## load the reference table
ref_oligo_med2 <- read.table("../../../../output/2_Oligotyping/2B/2B_REF_info_elizabethkingia.tsv", sep="\t", header = TRUE)

## select only the 2 oligotypes of Elizabethkingia
ref_oligo_med2 <- ref_oligo_med2[!is.na(ref_oligo_med2$oligotype),]

## change order of columns
ref_oligo_med2 <- ref_oligo_med2 %>% select(c(seq, oligotype, MED_node_frequency_size, OLIGO_oligotype_frequency_size))

## create a column with reference name (will be used in plots)
ref_oligo_med2$ref <- paste0("oligotype_", ref_oligo_med2$OLIGO_oligotype_frequency_size, " / node_", ref_oligo_med2$MED_node_frequency_size)

## create a copy of fasta 
fasta2 <- fasta

# Matrix count

## create an oligotype column in the matrix count
matrix_count$oligotype <- rownames(matrix_count)

## change order of columns
matrix_count <- matrix_count %>% dplyr::select(c(oligotype, everything()))

## merge the matrix count and the reference dataframe
matrix_count2 <- matrix_count %>% merge(ref_oligo_med2 %>% dplyr::select(-c(seq)), by="oligotype")

## change order of columns
matrix_count2 <- matrix_count2 %>% dplyr::select(c(oligotype, MED_node_frequency_size, OLIGO_oligotype_frequency_size, ref, everything()))

## change rownames
rownames(matrix_count2) <- matrix_count2$ref

## change order of columns
matrix_count2 <- matrix_count2 %>% dplyr::select(-c(oligotype, ref, MED_node_frequency_size, OLIGO_oligotype_frequency_size))

## print it
matrix_count2 %>%
  kbl() %>%
  kable_paper("hover", full_width = F)

	CTC1	CTC10	CTC11	CTC12	CTC13	CTC14	CTC15	CTC2	CTC3	CTC4	CTC5	CTC6	CTC7	CTC9	NP27	NP30	NP36	S126	S146	S147	S154	S160	S163	S164	S165	S166	S18	S19	S20	S21	S23	S24	S25	S26	S27	S28	S29	S30	S31	S32	S33	S34	S35	S36	S37	S38	S39	S40	S41	S44	S45	S46	S47	S48	S49	S50	S51	S52
oligotype_C (58) \| size:124095 / node_N1160 (58) \| size:115121	19626	592	5363	287	2687	1000	14770	1700	2520	1075	326	10704	15	2688	11	1	6	1	4	1	1	1	1	6	2	2	1	12104	12	6	33	53	3	8	1755	5	4	1933	4229	2020	2060	225	392	17944	3328	4744	4931	4388	1	36	32	9	6	221	43	71	61	47
oligotype_G (43) \| size:11092 / node_N0990 (43) \| size:11092	3017	80	963	43	437	150	2088	254	378	165	60	1649	5	399	1	0	1	0	0	0	0	0	0	1	0	4	0	277	0	0	5	1	0	1	27	0	1	55	121	36	66	5	10	390	71	73	110	130	0	1	1	0	1	6	2	3	1	3

## edit the fasta dataframe
fasta2 <- fasta2 %>% merge(ref_oligo_med2 %>% dplyr::select(-c(seq)), by="oligotype")
rownames(fasta2) <- fasta2$ref
fasta2 <- fasta2 %>% dplyr::select(-c(MED_node_frequency_size, OLIGO_oligotype_frequency_size, oligotype))

## print it
fasta2 %>%
  kbl() %>%
  kable_paper("hover", full_width = F)

	ref
oligotype_C (58) \| size:124095 / node_N1160 (58) \| size:115121	oligotype_C (58) \| size:124095 / node_N1160 (58) \| size:115121
oligotype_G (43) \| size:11092 / node_N0990 (43) \| size:11092	oligotype_G (43) \| size:11092 / node_N0990 (43) \| size:11092

Metadata

metadata <- read.csv("../../../../metadata/metadata.csv", sep=";")
rownames(metadata) <- metadata$Sample

Phyloseq object with oligotypes

# convert matrix_count into matrix and numeric
matrix_count <- matrix_count2 %>% as.matrix()
class(matrix_count) <- "numeric"

# phyloseq elements
OTU = otu_table(as.matrix(matrix_count), taxa_are_rows =TRUE)
TAX = tax_table(as.matrix(fasta2))
SAM = sample_data(metadata)

# phyloseq object
ps <- phyloseq(OTU, TAX, SAM)
ps

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 2 taxa and 58 samples ]
## sample_data() Sample Data:       [ 58 samples by 8 sample variables ]
## tax_table()   Taxonomy Table:    [ 2 taxa by 1 taxonomic ranks ]

compute_read_counts(ps)

## [1] 135187

# remove blanks
ps <- subset_samples(ps, Strain!="Blank")
ps <- check_ps(ps)
ps

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 2 taxa and 56 samples ]
## sample_data() Sample Data:       [ 56 samples by 8 sample variables ]
## tax_table()   Taxonomy Table:    [ 2 taxa by 1 taxonomic ranks ]

Create new metadata with Percent

Load ps with all samples (for final plot)

ps.filter <- readRDS("../../../../output/1_MED/1D/1D_MED_phyloseq_decontam.rds")
ps.filter <- check_ps(ps.filter)

Edit new metadata with Percent_elizabethkingia

guide_italics <- guides(fill = guide_legend(label.theme = element_text(size = 16, face = "italic", colour = "Black", angle = 0)))

# add read depth in sample table of phyloseq object
sample_data(ps.filter)$Read_depth <- sample_sums(ps.filter)

# select Elizabethkingia
ps.elizabethkingia <- ps.filter %>% subset_taxa(Genus=="Elizabethkingia")

# add read depth of Elizabethkingia
sample_data(ps.filter)$Read_elizabethkingia <- sample_sums(ps.elizabethkingia)
sample_data(ps.filter) %>% colnames()

##  [1] "Sample"               "Strain"               "Field"               
##  [4] "Country"              "Organ"                "Species"             
##  [7] "Run"                  "Control"              "Species_italic"      
## [10] "Strain_italic"        "Read_depth"           "is.neg"              
## [13] "Read_elizabethkingia"

sample_data(ps.elizabethkingia) %>% colnames()

##  [1] "Sample"         "Strain"         "Field"          "Country"       
##  [5] "Organ"          "Species"        "Run"            "Control"       
##  [9] "Species_italic" "Strain_italic"  "Read_depth"     "is.neg"

# add percent of Elizabethkingia
sample_data(ps.filter)$Percent_elizabethkingia <- sample_data(ps.filter)$Read_elizabethkingia / sample_data(ps.filter)$Read_depth

# round the percent of Elizabethkingia at 2 decimals
sample_data(ps.filter)$Percent_elizabethkingia <- sample_data(ps.filter)$Percent_elizabethkingia %>% round(2)

# extract metadata table
test <- data.frame(sample_data(ps.filter))

# merge this metadata table with the other
new.metadata <- data.frame(sample_data(ps)) %>% merge(test %>% dplyr::select(c(Sample, Read_depth, Read_elizabethkingia, Percent_elizabethkingia)), by="Sample")
new.metadata <- test[new.metadata$Sample %in% sample_names(ps),]
rownames(new.metadata) <- new.metadata$Sample

# print it
new.metadata %>%
  kbl() %>%
  kable_paper("hover", full_width = F)

	Sample	Strain	Field	Country	Organ	Species	Run	Control	Species_italic	Strain_italic	Read_depth	is.neg	Read_elizabethkingia	Percent_elizabethkingia
CTC1	CTC1	Laboratory - Slab TC (Wolbachia -)	Lab	France	Whole	Culex quinquefasciatus	run2	True sample	italic(“Culex quinquefasciatus”)	paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”)	29779	FALSE	22643	0.76
CTC10	CTC10	Laboratory - Slab TC (Wolbachia -)	Lab	France	Whole	Culex quinquefasciatus	run2	True sample	italic(“Culex quinquefasciatus”)	paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”)	2609	FALSE	672	0.26
CTC11	CTC11	Laboratory - Slab TC (Wolbachia -)	Lab	France	Whole	Culex quinquefasciatus	run2	True sample	italic(“Culex quinquefasciatus”)	paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”)	13874	FALSE	6326	0.46
CTC12	CTC12	Laboratory - Slab TC (Wolbachia -)	Lab	France	Whole	Culex quinquefasciatus	run2	True sample	italic(“Culex quinquefasciatus”)	paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”)	1146	FALSE	330	0.29
CTC13	CTC13	Laboratory - Slab TC (Wolbachia -)	Lab	France	Whole	Culex quinquefasciatus	run2	True sample	italic(“Culex quinquefasciatus”)	paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”)	18035	FALSE	3124	0.17
CTC14	CTC14	Laboratory - Slab TC (Wolbachia -)	Lab	France	Whole	Culex quinquefasciatus	run2	True sample	italic(“Culex quinquefasciatus”)	paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”)	1708	FALSE	1150	0.67
CTC15	CTC15	Laboratory - Slab TC (Wolbachia -)	Lab	France	Whole	Culex quinquefasciatus	run2	True sample	italic(“Culex quinquefasciatus”)	paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”)	23180	FALSE	16858	0.73
CTC2	CTC2	Laboratory - Slab TC (Wolbachia -)	Lab	France	Whole	Culex quinquefasciatus	run2	True sample	italic(“Culex quinquefasciatus”)	paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”)	30692	FALSE	1954	0.06
CTC3	CTC3	Laboratory - Slab TC (Wolbachia -)	Lab	France	Whole	Culex quinquefasciatus	run2	True sample	italic(“Culex quinquefasciatus”)	paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”)	39920	FALSE	2898	0.07
CTC4	CTC4	Laboratory - Slab TC (Wolbachia -)	Lab	France	Whole	Culex quinquefasciatus	run2	True sample	italic(“Culex quinquefasciatus”)	paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”)	2139	FALSE	1240	0.58
CTC5	CTC5	Laboratory - Slab TC (Wolbachia -)	Lab	France	Whole	Culex quinquefasciatus	run2	True sample	italic(“Culex quinquefasciatus”)	paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”)	15789	FALSE	386	0.02
CTC6	CTC6	Laboratory - Slab TC (Wolbachia -)	Lab	France	Whole	Culex quinquefasciatus	run2	True sample	italic(“Culex quinquefasciatus”)	paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”)	19753	FALSE	12353	0.63
CTC9	CTC9	Laboratory - Slab TC (Wolbachia -)	Lab	France	Whole	Culex quinquefasciatus	run2	True sample	italic(“Culex quinquefasciatus”)	paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”)	4980	FALSE	3087	0.62
NP14	NP14	Field - Guadeloupe	Field	Guadeloupe	Ovary	Aedes aegypti	run3	True sample	italic(“Aedes aegypti”)	Field_-Guadeloupe	7973	FALSE	0	0.00
NP2	NP2	Field - Guadeloupe	Field	Guadeloupe	Ovary	Culex quinquefasciatus	run3	True sample	italic(“Culex quinquefasciatus”)	Field_-Guadeloupe	648335	FALSE	0	0.00
NP20	NP20	Field - Guadeloupe	Field	Guadeloupe	Ovary	Aedes aegypti	run3	True sample	italic(“Aedes aegypti”)	Field_-Guadeloupe	136	FALSE	0	0.00
NP27	NP27	Field - Guadeloupe	Field	Guadeloupe	Whole	Culex quinquefasciatus	run3	True sample	italic(“Culex quinquefasciatus”)	Field_-Guadeloupe	1234	FALSE	12	0.01
NP29	NP29	Field - Guadeloupe	Field	Guadeloupe	Whole	Culex quinquefasciatus	run3	True sample	italic(“Culex quinquefasciatus”)	Field_-Guadeloupe	203	FALSE	0	0.00
NP30	NP30	Field - Guadeloupe	Field	Guadeloupe	Whole	Culex quinquefasciatus	run3	True sample	italic(“Culex quinquefasciatus”)	Field_-Guadeloupe	228	FALSE	1	0.00
NP34	NP34	Field - Guadeloupe	Field	Guadeloupe	Whole	Culex quinquefasciatus	run3	True sample	italic(“Culex quinquefasciatus”)	Field_-Guadeloupe	95	FALSE	0	0.00
NP35	NP35	Field - Guadeloupe	Field	Guadeloupe	Whole	Aedes aegypti	run3	True sample	italic(“Aedes aegypti”)	Field_-Guadeloupe	196532	FALSE	0	0.00
NP36	NP36	Field - Guadeloupe	Field	Guadeloupe	Whole	Aedes aegypti	run3	True sample	italic(“Aedes aegypti”)	Field_-Guadeloupe	249	FALSE	7	0.03
NP37	NP37	Field - Guadeloupe	Field	Guadeloupe	Whole	Aedes aegypti	run3	True sample	italic(“Aedes aegypti”)	Field_-Guadeloupe	419340	FALSE	0	0.00
NP38	NP38	Field - Guadeloupe	Field	Guadeloupe	Whole	Aedes aegypti	run3	True sample	italic(“Aedes aegypti”)	Field_-Guadeloupe	282479	FALSE	0	0.00
NP39	NP39	Field - Guadeloupe	Field	Guadeloupe	Whole	Aedes aegypti	run3	True sample	italic(“Aedes aegypti”)	Field_-Guadeloupe	218684	FALSE	0	0.00
NP41	NP41	Field - Guadeloupe	Field	Guadeloupe	Whole	Aedes aegypti	run3	True sample	italic(“Aedes aegypti”)	Field_-Guadeloupe	247152	FALSE	0	0.00
NP42	NP42	Field - Guadeloupe	Field	Guadeloupe	Whole	Aedes aegypti	run3	True sample	italic(“Aedes aegypti”)	Field_-Guadeloupe	185157	FALSE	0	0.00
NP43	NP43	Field - Guadeloupe	Field	Guadeloupe	Whole	Aedes aegypti	run3	True sample	italic(“Aedes aegypti”)	Field_-Guadeloupe	239335	FALSE	0	0.00
NP44	NP44	Field - Guadeloupe	Field	Guadeloupe	Whole	Aedes aegypti	run3	True sample	italic(“Aedes aegypti”)	Field_-Guadeloupe	156879	FALSE	0	0.00
NP5	NP5	Field - Guadeloupe	Field	Guadeloupe	Ovary	Culex quinquefasciatus	run3	True sample	italic(“Culex quinquefasciatus”)	Field_-Guadeloupe	736159	FALSE	0	0.00
NP8	NP8	Field - Guadeloupe	Field	Guadeloupe	Ovary	Culex quinquefasciatus	run3	True sample	italic(“Culex quinquefasciatus”)	Field_-Guadeloupe	334799	FALSE	0	0.00
S100	S100	Field - Camping Europe	Field	France	Ovary	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Camping~Europe	52486	FALSE	0	0.00
S102	S102	Field - Camping Europe	Field	France	Ovary	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Camping~Europe	3456	FALSE	0	0.00
S104	S104	Field - Camping Europe	Field	France	Ovary	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Camping~Europe	52403	FALSE	0	0.00
S105	S105	Field - Camping Europe	Field	France	Ovary	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Camping~Europe	55577	FALSE	0	0.00
S106	S106	Field - Camping Europe	Field	France	Ovary	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Camping~Europe	33053	FALSE	0	0.00
S107	S107	Field - Camping Europe	Field	France	Ovary	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Camping~Europe	52154	FALSE	0	0.00
S108	S108	Field - Camping Europe	Field	France	Ovary	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Camping~Europe	55735	FALSE	0	0.00
S109	S109	Field - Camping Europe	Field	France	Ovary	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Camping~Europe	59023	FALSE	0	0.00
S110	S110	Field - Camping Europe	Field	France	Ovary	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Camping~Europe	57377	FALSE	0	0.00
S121	S121	Field - Bosc	Field	France	Ovary	Culex pipiens	run2	True sample	italic(“Culex pipiens”)	Field_-Bosc	20361	FALSE	0	0.00
S122	S122	Field - Bosc	Field	France	Ovary	Culex pipiens	run2	True sample	italic(“Culex pipiens”)	Field_-Bosc	9803	FALSE	0	0.00
S123	S123	Field - Bosc	Field	France	Ovary	Culex pipiens	run2	True sample	italic(“Culex pipiens”)	Field_-Bosc	20130	FALSE	0	0.00
S124	S124	Field - Bosc	Field	France	Ovary	Culex pipiens	run2	True sample	italic(“Culex pipiens”)	Field_-Bosc	18146	FALSE	0	0.00
S126	S126	Field - Bosc	Field	France	Ovary	Culex pipiens	run2	True sample	italic(“Culex pipiens”)	Field_-Bosc	15235	FALSE	1	0.00
S127	S127	Field - Bosc	Field	France	Ovary	Culex pipiens	run2	True sample	italic(“Culex pipiens”)	Field_-Bosc	24696	FALSE	0	0.00
S128	S128	Field - Bosc	Field	France	Ovary	Culex pipiens	run2	True sample	italic(“Culex pipiens”)	Field_-Bosc	16305	FALSE	0	0.00
S146	S146	Laboratory - Lavar	Lab	France	Ovary	Culex pipiens	run2	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	25012	FALSE	4	0.00
S147	S147	Laboratory - Lavar	Lab	France	Ovary	Culex pipiens	run2	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	25171	FALSE	1	0.00
S148	S148	Laboratory - Lavar	Lab	France	Ovary	Culex pipiens	run2	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	14164	FALSE	0	0.00
S150	S150	Laboratory - Lavar	Lab	France	Ovary	Culex pipiens	run2	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	15081	FALSE	0	0.00
S151	S151	Laboratory - Lavar	Lab	France	Ovary	Culex pipiens	run2	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	22944	FALSE	0	0.00
S152	S152	Laboratory - Lavar	Lab	France	Ovary	Culex pipiens	run2	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	15082	FALSE	0	0.00
S153	S153	Laboratory - Lavar	Lab	France	Ovary	Culex pipiens	run2	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	17040	FALSE	0	0.00
S154	S154	Laboratory - Lavar	Lab	France	Ovary	Culex pipiens	run2	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	9626	FALSE	1	0.00
S160	S160	Laboratory - Lavar	Lab	France	Ovary	Culex pipiens	run2	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	72508	FALSE	1	0.00
S162	S162	Laboratory - Lavar	Lab	France	Ovary	Culex pipiens	run2	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	25180	FALSE	0	0.00
S163	S163	Laboratory - Lavar	Lab	France	Ovary	Culex pipiens	run2	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	12333	FALSE	1	0.00
S164	S164	Laboratory - Lavar	Lab	France	Ovary	Culex pipiens	run2	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	22368	FALSE	7	0.00
S165	S165	Laboratory - Lavar	Lab	France	Ovary	Culex pipiens	run2	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	17731	FALSE	2	0.00
S166	S166	Field - Camping Europe	Field	France	Ovary	Culex pipiens	run2	True sample	italic(“Culex pipiens”)	Field_-Camping~Europe	13979	FALSE	6	0.00
S167	S167	Field - Bosc	Field	France	Ovary	Culex pipiens	run2	True sample	italic(“Culex pipiens”)	Field_-Bosc	14048	FALSE	0	0.00
S169	S169	Field - Camping Europe	Field	France	Ovary	Culex pipiens	run2	True sample	italic(“Culex pipiens”)	Field_-Camping~Europe	11553	FALSE	0	0.00
S170	S170	Field - Camping Europe	Field	France	Ovary	Culex pipiens	run2	True sample	italic(“Culex pipiens”)	Field_-Camping~Europe	8852	FALSE	0	0.00
S18	S18	Laboratory - Lavar	Lab	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	4290	FALSE	1	0.00
S19	S19	Laboratory - Lavar	Lab	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	44527	FALSE	12381	0.28
S20	S20	Laboratory - Lavar	Lab	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	42864	FALSE	12	0.00
S21	S21	Laboratory - Lavar	Lab	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	33798	FALSE	6	0.00
S22	S22	Laboratory - Lavar	Lab	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	19044	FALSE	0	0.00
S23	S23	Laboratory - Lavar	Lab	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	38172	FALSE	38	0.00
S24	S24	Laboratory - Lavar	Lab	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	42355	FALSE	54	0.00
S25	S25	Laboratory - Lavar	Lab	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	47688	FALSE	3	0.00
S26	S26	Laboratory - Lavar	Lab	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	5394	FALSE	9	0.00
S27	S27	Laboratory - Lavar	Lab	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	24558	FALSE	1782	0.07
S28	S28	Laboratory - Lavar	Lab	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	4503	FALSE	5	0.00
S30	S30	Laboratory - Lavar	Lab	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	25353	FALSE	1988	0.08
S31	S31	Laboratory - Lavar	Lab	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	20417	FALSE	4350	0.21
S32	S32	Laboratory - Lavar	Lab	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	12441	FALSE	2056	0.17
S33	S33	Laboratory - Lavar	Lab	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	33867	FALSE	2126	0.06
S34	S34	Laboratory - Lavar	Lab	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	9367	FALSE	230	0.02
S35	S35	Laboratory - Lavar	Lab	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	11663	FALSE	402	0.03
S36	S36	Laboratory - Lavar	Lab	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	33020	FALSE	18334	0.56
S37	S37	Laboratory - Lavar	Lab	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	18340	FALSE	3399	0.19
S38	S38	Laboratory - Lavar	Lab	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	54790	FALSE	4817	0.09
S39	S39	Laboratory - Lavar	Lab	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	36273	FALSE	5041	0.14
S40	S40	Laboratory - Lavar	Lab	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Laboratory_-Lavar	44448	FALSE	4518	0.10
S42	S42	Field - Camping Europe	Field	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Camping~Europe	4107	FALSE	0	0.00
S43	S43	Field - Camping Europe	Field	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Camping~Europe	9279	FALSE	0	0.00
S44	S44	Field - Camping Europe	Field	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Camping~Europe	8026	FALSE	37	0.00
S45	S45	Field - Camping Europe	Field	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Camping~Europe	18150	FALSE	33	0.00
S47	S47	Field - Camping Europe	Field	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Camping~Europe	1951	FALSE	7	0.00
S48	S48	Field - Camping Europe	Field	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Camping~Europe	56738	FALSE	227	0.00
S49	S49	Field - Bosc	Field	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Bosc	33498	FALSE	45	0.00
S50	S50	Field - Bosc	Field	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Bosc	28481	FALSE	74	0.00
S51	S51	Field - Bosc	Field	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Bosc	61788	FALSE	62	0.00
S52	S52	Field - Bosc	Field	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Bosc	21553	FALSE	50	0.00
S55	S55	Field - Bosc	Field	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Bosc	50447	FALSE	0	0.00
S56	S56	Field - Bosc	Field	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Bosc	42609	FALSE	0	0.00
S57	S57	Field - Bosc	Field	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Bosc	49157	FALSE	0	0.00
S58	S58	Field - Bosc	Field	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Bosc	30357	FALSE	0	0.00
S59	S59	Field - Bosc	Field	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Bosc	32798	FALSE	0	0.00
S60	S60	Field - Bosc	Field	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Bosc	44485	FALSE	0	0.00
S61	S61	Field - Bosc	Field	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Bosc	49545	FALSE	0	0.00
S63	S63	Field - Bosc	Field	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Bosc	53444	FALSE	0	0.00
S64	S64	Field - Bosc	Field	France	Whole	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Bosc	47628	FALSE	0	0.00
S79	S79	Field - Camping Europe	Field	France	Ovary	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Camping~Europe	59755	FALSE	0	0.00
S80	S80	Field - Camping Europe	Field	France	Ovary	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Camping~Europe	52788	FALSE	0	0.00
S83	S83	Field - Camping Europe	Field	France	Ovary	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Camping~Europe	42272	FALSE	0	0.00
S84	S84	Field - Camping Europe	Field	France	Ovary	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Camping~Europe	56676	FALSE	0	0.00
S85	S85	Field - Camping Europe	Field	France	Ovary	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Camping~Europe	41690	FALSE	0	0.00
S86	S86	Field - Camping Europe	Field	France	Ovary	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Camping~Europe	61984	FALSE	0	0.00
S87	S87	Field - Bosc	Field	France	Ovary	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Bosc	65958	FALSE	0	0.00
S88	S88	Field - Bosc	Field	France	Ovary	Culex pipiens	run1	True sample	italic(“Culex pipiens”)	Field_-Bosc	53102	FALSE	0	0.00

# replace metadata in the created phyloseq object
sample_data(ps) <- sample_data(new.metadata)

Taxonomic structure

Count

col <- brewer.pal(7, "Pastel2")

# reshape data for plot
test3 <- test %>% select(c(Sample, Species, Strain, Organ, Read_depth, Read_elizabethkingia)) %>% reshape2::melt(id.vars=c("Sample", "Species", "Strain", "Organ"), vars=c("Read_depth", "Read_elizabethkingia"))

count_whole <- test3[test3$Organ=="Whole",]
count_ovary <- test3[test3$Organ=="Ovary",]

make.italic <- function(x) as.expression(lapply(x, function(y) bquote(italic(.(y)))))

levels(count_whole$Species)= c("Aedes aegypti"=make.italic("Aedes aegypti"),
               "Culex pipiens"=make.italic("Culex pipiens"),
               "Culex quinquefasciatus"=make.italic("Culex quinquefasciatus"))

levels(count_ovary$Species)= c("Aedes aegypti"=make.italic("Aedes aegypti"),
               "Culex pipiens"=make.italic("Culex pipiens"),
               "Culex quinquefasciatus"=make.italic("Culex quinquefasciatus"))

levels(count_whole$Strain) <- c("Bosc", "Camping~Europe", "Guadeloupe", "Lavar~(lab)", expression(paste(italic("Wolbachia"), "- (Slab TC)")))

levels(count_ovary$Strain) <- c("Bosc", "Camping~Europe", "Guadeloupe", "Lavar~(lab)", expression(paste(italic("Wolbachia"), "- (Slab TC)")))


# plot
p_count1 <- ggplot(count_whole, aes(x = Sample, y = value, fill=variable))+ 
  geom_bar(position = "dodge", stat = "identity")+
  scale_fill_manual(values = col)+
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, size=12, hjust=1, vjust=0.5)) +
  ggtitle("") + 
  guide_italics+
  theme(legend.title = element_text(size = 20), 
        legend.position="bottom",
        legend.text=element_text(size=14), 
        panel.spacing.y=unit(1, "lines"), 
        panel.spacing.x=unit(0.8, "lines"),
        panel.spacing=unit(0,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        strip.text.x = element_text(size = 16),
        panel.border=element_rect(color="grey90"),
        axis.ticks.x=element_blank(),
        axis.text.y = element_text(size=18)) +
  facet_wrap(~Species+Strain+Organ, scales = "free_x", ncol=3, labeller=label_parsed)+
  labs(y="Sequence counts")+
  ylim(0, 900000)+
  geom_text(aes(label=value), position=position_dodge(width=1.1), width=0.25, size=4, hjust=-0.25, vjust=0.5, angle=90)+
  guides(fill=guide_legend(title="Read"))

## Warning: Ignoring unknown parameters: width

p_count2 <- ggplot(count_ovary, aes(x = Sample, y = value, fill=variable))+ 
  geom_bar(position = "dodge", stat = "identity")+
    scale_fill_manual(values = col)+
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, size=18, hjust=1, vjust=0.5)) +
  ggtitle("") + 
  guide_italics+
  theme(legend.title = element_text(size = 20), 
        legend.position="bottom",
        legend.text=element_text(size=14), 
        panel.spacing.y=unit(1, "lines"), 
        panel.spacing.x=unit(0.8, "lines"),
        panel.spacing=unit(0,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        strip.text.x = element_text(size = 16),
        panel.border=element_rect(color="grey90"),
        axis.ticks.x=element_blank(),
        axis.text.y = element_text(size=18)) +
 facet_wrap(~Species+Strain+Organ, scales = "free_x", ncol=3, labeller=label_parsed)+
  labs(y="Sequence counts")+
    ylim(0, 900000)+
  geom_text(aes(label=value), position=position_dodge(width=0.8), width=0.25, size=4, hjust=-0.25, vjust=0.5, angle=90)+
  guides(fill=guide_legend(title="Read"))

## Warning: Ignoring unknown parameters: width

# afficher plot
p_count1

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

p_count2

# panels
p_group <- plot_grid(p_count1+theme(legend.position="none"), 
          p_count2+theme(legend.position="none"), 
          nrow=2, 
          ncol=1)+
    draw_plot_label(c("B1", "B2"), c(0, 0), c(1, 0.5), size = 20)

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

legend_plot <- get_legend(p_count1 + theme(legend.position="bottom"))

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

p_counts <- plot_grid(p_group, legend_plot, nrow=2, ncol=1, rel_heights = c(1, .1))
p_counts

Whole (the most abundant nodes)

guide_italics <- guides(fill = guide_legend(label.theme = element_text(size = 16, face = "italic", colour = "Black", angle = 0),
                                            nrow=2, byrow=TRUE))

# select whole
ps.filter.whole <- subset_samples(ps, Organ=="Whole")
ps.filter.whole <- prune_taxa(taxa_sums(ps.filter.whole) >= 1, ps.filter.whole)
ps.filter.whole <- prune_samples(sample_sums(ps.filter.whole) >= 1, ps.filter.whole)
ps.filter.whole

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 2 taxa and 45 samples ]
## sample_data() Sample Data:       [ 45 samples by 14 sample variables ]
## tax_table()   Taxonomy Table:    [ 2 taxa by 1 taxonomic ranks ]

# data pour plot
#data_for_plot2 <- taxo_data_fast(ps.filter.whole, method = "abundance")
data_for_plot2 <- taxo_data(ps.filter.whole, top=15)

## Warning in psmelt(ps_global): The sample variables: 
## Sample
##  have been renamed to: 
## sample_Sample
## to avoid conflicts with special phyloseq plot attribute names.

paste0("\n15 MOST ABUNDANT GENUS: \n") %>% cat()

## 
## 15 MOST ABUNDANT GENUS:

paste0("\"", levels(data_for_plot2$Name), "\",\n") %>% cat()

## "oligotype_C (58) | size:124095 / node_N1160 (58) | size:115121.",
##  "oligotype_G (43) | size:11092 / node_N0990 (43) | size:11092.",

data_for_plot2$Name <- data_for_plot2$Name %>% gsub(pattern = "node_", replacement ="" ) %>% as.factor()
data_for_plot2$Name <- as.factor(data_for_plot2$Name)

new_names <- c("oligotype_C (58) | size:124095 / N1160 (58) | size:115121.",
               "oligotype_G (43) | size:11092 / N0990 (43) | size:11092.")

data_for_plot2$Name <- factor(data_for_plot2$Name, levels = new_names)

col_add <- brewer.pal(8, "Accent")

col <- c("oligotype_C (58) | size:124095 / N1160 (58) | size:115121."="#6AEEF7",
         "oligotype_G (43) | size:11092 / N0990 (43) | size:11092."="#6AC3F7")

levels(data_for_plot2$Species)= c("Aedes aegypti"=make.italic("Aedes aegypti"),
               "Culex pipiens"=make.italic("Culex pipiens"),
               "Culex quinquefasciatus"=make.italic("Culex quinquefasciatus"))

levels(data_for_plot2$Strain) <- c("Bosc", "Camping~Europe", "Guadeloupe", "Lavar~(lab)", expression(paste(italic("Wolbachia"), "- (Slab TC)")))

#data_for_plot2 <- data_for_plot2 %>% na.omit()

p2 <- ggplot(data_for_plot2, aes(x = Sample, y = Relative_Abundance, fill = Name, species=Species, organ=Organ, Strain=Strain))+ 
  geom_bar(position = "stack", stat = "identity")+
  scale_fill_manual(values = col)+
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, size=18, hjust=1, vjust=0.5)) +
  ggtitle("") + 
  guide_italics+
  theme(legend.title = element_text(size = 20),
        legend.position="bottom",
        legend.text = element_text(size=14),
        #legend.key.height = unit(1, 'cm'),
        panel.spacing.y=unit(1, "lines"),
        panel.spacing.x=unit(0.8, "lines"),
        panel.spacing=unit(0,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        strip.text.x = element_text(size = 16),
        panel.border=element_rect(color="grey90"),
        axis.ticks.x=element_blank(),
        axis.text.y = element_text(size=18)) +
  facet_wrap(~Species+Strain+Organ, scales = "free", ncol=3, labeller=label_parsed)+
  labs(x="Sample", y="Relative abundance", fill="Oligotype / MED node")

p2

Ovary (the most abundant nodes)

# select ovary
ps.filter.ovary <- subset_samples(ps, Organ=="Ovary")
ps.filter.ovary <- prune_taxa(taxa_sums(ps.filter.ovary) >= 1, ps.filter.ovary)
ps.filter.ovary <- prune_samples(sample_sums(ps.filter.ovary) >= 1, ps.filter.ovary)
ps.filter.ovary

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 2 taxa and 9 samples ]
## sample_data() Sample Data:       [ 9 samples by 14 sample variables ]
## tax_table()   Taxonomy Table:    [ 2 taxa by 1 taxonomic ranks ]

# data pour plot
#data_for_plot3 <- taxo_data_fast(ps.filter.ovary, method = "abundance")
data_for_plot3 <- taxo_data(ps.filter.ovary, top=15)

## Warning in psmelt(ps_global): The sample variables: 
## Sample
##  have been renamed to: 
## sample_Sample
## to avoid conflicts with special phyloseq plot attribute names.

paste0("\n15 MOST ABUNDANT GENUS: \n") %>% cat()

## 
## 15 MOST ABUNDANT GENUS:

paste0("\"", levels(data_for_plot3$Name), "\",\n") %>% cat()

## "oligotype_C (58) | size:124095 / node_N1160 (58) | size:115121.",
##  "oligotype_G (43) | size:11092 / node_N0990 (43) | size:11092.",

data_for_plot3$Name <- data_for_plot3$Name %>% gsub(pattern = "node_", replacement ="" ) %>% as.factor()
data_for_plot3$Name <- as.factor(data_for_plot3$Name)

new_names <- c("oligotype_C (58) | size:124095 / N1160 (58) | size:115121.",
               "oligotype_G (43) | size:11092 / N0990 (43) | size:11092.")

data_for_plot3$Name <- factor(data_for_plot3$Name, levels = new_names)

# levels(data_for_plot3$Species)= c("Aedes aegypti"=make.italic("Aedes aegypti"),
#                "Culex pipiens"=make.italic("Culex pipiens"),
#                "Culex quinquefasciatus"=make.italic("Culex quinquefasciatus"))

levels(data_for_plot3$Species)= c("Culex pipiens"=make.italic("Culex pipiens"))

#levels(data_for_plot3$Strain) <- c("Bosc", "Camping~Europe", "Guadeloupe", "Lavar~(lab)", expression(paste(italic("Wolbachia"), "- (Slab TC)")))
levels(data_for_plot3$Strain) <- c("Bosc", "Camping~Europe", "Lavar~(lab)")

#data_for_plot3 <- data_for_plot3 %>% na.omit()

p3 <- ggplot(data_for_plot3, aes(x = Sample, y = Relative_Abundance, fill = Name, species=Species, organ=Organ, Strain=Strain))+ 
  geom_bar(position = "stack", stat = "identity")+
  scale_fill_manual(values = col)+
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, size=18, hjust=1, vjust=0.5)) +
  ggtitle("") + 
  guide_italics+
  theme(legend.title = element_text(size = 20), 
        legend.position="bottom",
        legend.text = element_text(size=14),
        #legend.key.height = unit(1, 'cm'),
        panel.spacing.y=unit(1, "lines"), 
        panel.spacing.x=unit(0.8, "lines"),
        panel.spacing=unit(0,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        strip.text.x = element_text(size = 16),
        panel.border=element_rect(color="grey90"),
        axis.ticks.x=element_blank(),
        axis.text.y = element_text(size=18)) +
  facet_wrap(~Species+Strain+Organ, scales = "free", ncol=3, labeller=label_parsed)+
  labs(x="Sample", y="Relative abundance", fill="Oligotype / MED node")

p3

Panels taxonomy of whole / ovary

legend_plot <- get_legend(p2 + theme(legend.position="bottom"))

# panels
p_group <- plot_grid(p2+theme(legend.position="none"), 
          p3+theme(legend.position="none"), 
          nrow=2, 
          ncol=1)+
    draw_plot_label(c("A1", "A2"), c(0, 0), c(1, 0.5), size = 20)

p_taxo <- plot_grid(p_group, legend_plot, nrow=2, rel_heights = c(1, .1))
p_taxo

Save taxonomic plot

tiff("../../../../output/2_Oligotyping/2D/2Dd_OLIGO_counts_elizabethkingia.tiff", units="in", width=20, height=18, res=300)
p_counts
dev.off()

## quartz_off_screen 
##                 2

tiff("../../../../output/2_Oligotyping/2D/2Dd_OLIGO_taxonomic_elizabethkingia_whole.tiff", units="in", width=16, height=12, res=300)
p2
dev.off()

## quartz_off_screen 
##                 2

tiff("../../../../output/2_Oligotyping/2D/2Dd_OLIGO_taxonomic_elizabethkingia_ovary.tiff", units="in", width=18, height=14, res=300)
p3
dev.off()

## quartz_off_screen 
##                 2

tiff("../../../../output/2_Oligotyping/2D/2Dd_OLIGO_taxonomic_elizabethkingia.tiff", units="in", width=18, height=16, res=300)
p_taxo
dev.off()

## quartz_off_screen 
##                 2

png("../../../../output/2_Oligotyping/2D/2Dd_OLIGO_counts_elizabethkingia_big.png", units="in", width=20, height=18, res=300)
p_counts
dev.off()

## quartz_off_screen 
##                 2

png("../../../../output/2_Oligotyping/2D/2Dd_OLIGO_counts_elizabethkingia_small.png", units="in", width=18, height=14, res=300)
p_counts
dev.off()

## quartz_off_screen 
##                 2

png("../../../../output/2_Oligotyping/2D/2Dd_OLIGO_taxonomic_elizabethkingia_whole.png", units="in", width=16, height=12, res=300)
p2
dev.off()

## quartz_off_screen 
##                 2

png("../../../../output/2_Oligotyping/2D/2Dd_OLIGO_taxonomic_elizabethkingia_ovary.png", units="in", width=18, height=14, res=300)
p3
dev.off()

## quartz_off_screen 
##                 2

png("../../../../output/2_Oligotyping/2D/2Dd_OLIGO_taxonomic_elizabethkingia_big.png", units="in", width=18, height=18, res=300)
p_taxo
dev.off()

## quartz_off_screen 
##                 2

png("../../../../output/2_Oligotyping/2D/2Dd_OLIGO_taxonomic_elizabethkingia_small.png", units="in", width=18, height=14, res=300)
p_taxo
dev.off()

## quartz_off_screen 
##                 2

Make main plot

img <- magick::image_read(paste0(path_elizabeth, "/HTML-OUTPUT/entropy.png"))
p_entropy <- magick::image_ggplot(img, interpolate = TRUE)
p_entropy+ theme(plot.margin = unit(c(-7,-2.5,-7,-0.5), "cm"))

p_entropy+ theme(plot.margin=unit(c(-7,-2,-12,-5), "mm"))

aligned <- plot_grid(p_taxo, 
                     p_counts, 
                     align="hv")

aligned

p_entropy2 <- plot_grid(p_entropy, nrow=1)+
  draw_plot_label(c("C"), c(0), c(1), size=20, hjust=-0.5)

p_entropy2

t_plot <- plot_grid(aligned, 
                    p_entropy2,
                    nrow=2, 
                    ncol=1, 
                    scale=1,
                    rel_heights=c(2,1))

t_plot

tiff("../../../../output/2_Oligotyping/2D/2Dd_OLIGO_main_elizabethkingia.tiff", width=36, height=36, res=300, units="in")
t_plot
dev.off()

## quartz_off_screen 
##                 2

png("../../../../output/2_Oligotyping/2D/2Dd_OLIGO_main_elizabethkingia.png", width=36, height=36, res=300, units="in")
t_plot
dev.off()

## quartz_off_screen 
##                 2

Session info

sessionInfo()

## R version 3.6.3 (2020-02-29)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS  10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] fr_FR.UTF-8/fr_FR.UTF-8/fr_FR.UTF-8/C/fr_FR.UTF-8/fr_FR.UTF-8
## 
## attached base packages:
## [1] stats4    parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] kableExtra_1.3.1    openxlsx_4.2.3      Biostrings_2.54.0  
##  [4] XVector_0.26.0      IRanges_2.20.2      S4Vectors_0.24.4   
##  [7] BiocGenerics_0.32.0 cowplot_1.1.0       ggpubr_0.4.0       
## [10] RColorBrewer_1.1-2  forcats_0.5.0       stringr_1.4.0      
## [13] dplyr_1.0.2         purrr_0.3.4         readr_1.4.0        
## [16] tidyr_1.1.2         tibble_3.0.4        tidyverse_1.3.0    
## [19] ggplot2_3.3.2       phyloseq_1.30.0    
## 
## loaded via a namespace (and not attached):
##  [1] colorspace_2.0-0  ggsignif_0.6.0    ellipsis_0.3.1    rio_0.5.16       
##  [5] fs_1.5.0          rstudioapi_0.11   farver_2.0.3      fansi_0.4.1      
##  [9] lubridate_1.7.9   xml2_1.3.2        codetools_0.2-16  splines_3.6.3    
## [13] knitr_1.30        ade4_1.7-15       jsonlite_1.7.1    broom_0.7.2      
## [17] cluster_2.1.0     dbplyr_1.4.4      compiler_3.6.3    httr_1.4.2       
## [21] backports_1.1.10  assertthat_0.2.1  Matrix_1.2-18     cli_2.1.0        
## [25] htmltools_0.5.1.1 tools_3.6.3       igraph_1.2.6      gtable_0.3.0     
## [29] glue_1.4.2        reshape2_1.4.4    Rcpp_1.0.5        carData_3.0-4    
## [33] Biobase_2.46.0    cellranger_1.1.0  jquerylib_0.1.3   vctrs_0.3.4      
## [37] multtest_2.42.0   ape_5.4-1         nlme_3.1-149      iterators_1.0.13 
## [41] xfun_0.22         ps_1.4.0          rvest_0.3.6       lifecycle_0.2.0  
## [45] rstatix_0.6.0     zlibbioc_1.32.0   MASS_7.3-53       scales_1.1.1     
## [49] hms_0.5.3         biomformat_1.14.0 rhdf5_2.30.1      yaml_2.2.1       
## [53] curl_4.3          sass_0.3.1        stringi_1.5.3     highr_0.8        
## [57] foreach_1.5.1     permute_0.9-5     zip_2.1.1         rlang_0.4.10     
## [61] pkgconfig_2.0.3   evaluate_0.14     lattice_0.20-41   Rhdf5lib_1.8.0   
## [65] labeling_0.4.2    tidyselect_1.1.0  plyr_1.8.6        magrittr_1.5     
## [69] bookdown_0.22     R6_2.4.1          magick_2.5.2      generics_0.0.2   
## [73] DBI_1.1.0         pillar_1.4.6      haven_2.3.1       foreign_0.8-75   
## [77] withr_2.3.0       mgcv_1.8-33       survival_3.2-7    abind_1.4-5      
## [81] modelr_0.1.8      crayon_1.3.4      car_3.0-10        rmarkdown_2.7    
## [85] grid_3.6.3        readxl_1.3.1      data.table_1.13.2 blob_1.2.1       
## [89] vegan_2.5-6       rmdformats_1.0.2  reprex_0.3.0      digest_0.6.26    
## [93] webshot_0.5.2     munsell_0.5.0     viridisLite_0.3.0 bslib_0.2.4